%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 遗传算法
% 简介：遗传算法寻找最大值
% 作者：Zhaojiang
% 日期：2023/10/9
% 企鹅：277746470
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;clear;clc;
%% target fucntion
[x,y]=meshgrid(-4:0.1:4);
z=target(x,y);
mesh(x,y,z);
title('Target Function')
xlabel('x');ylabel('y');zlabel('z');
pause(2);
close all;clear;clc;
%% 设定算法参数
rng(1234);
var_num=2;% 2个求解变量
target_resolution=1e-4;% 解的分辨率（小数后4位）
upper_bound = 4;% 解的最大值（小数后4位）
lower_bound = -4;% 解的最小值（小数后4位）
elitism = true;% 精英选择
population_size = 200;% 种群大小
generation_size = 100;% 迭代次数
cross_rate = 0.3;% 交叉概率
mutate_rate = 0.01;% 变异概率
% 计算染色体长度（二进制编码）
chromosome_size = 1;
bound_length = upper_bound-lower_bound;
while(1)
    resolution = bound_length/(2^(chromosome_size-1)-1);
    next_bit_resolution = bound_length/(2^chromosome_size-1);
    % 前一个长度不符合精度，这个长度符合精度
    if  (resolution >= target_resolution)&& ...
        (next_bit_resolution <= target_resolution)
        break;
    end
    chromosome_size=chromosome_size+1;
end
var_chromosome_size = chromosome_size;% 每个变量染色体长度
chromosome_size = chromosome_size * var_num;
%% 算法数据初始化
var_value = zeros(1,var_num);% 变量值(十进制)
fitness_value = zeros(1,population_size);% 适应度值
fitness_avg = zeros(1,generation_size);% 平均适应度
fitness_sum = zeros(1,population_size);% 总适应度
fitness_best = 0;% 当前最佳适应度
best_fitness = zeros(1,generation_size);% 记录最佳适应度
best_generation = 0;% 最佳繁殖代数
population =round(rand(population_size,chromosome_size));% 随机初始化种群
%% 开始迭代求解
for generation=1:generation_size
    % 计算适应度
    % 变量归零
    var_value = zeros(1,var_num);% 变量值(十进制)
    for i = 1:population_size
        % 种群染色体解码
        for j = 1:var_num
            for k = 1:var_chromosome_size 
                if population(i,k+(j-1)*var_chromosome_size) ~=0
                    var_value(j) = var_value(j) + 2^(k-1); 
                end
            end
            var_value(j) = lower_bound + var_value(j)*...
                        bound_length/(2^var_chromosome_size-1);
        end
        % 适应度计算（根据实际情况调整）
        fitness_value(i) = target(var_value(1),var_value(2));
    end
    % 种群排序
%     % 种群总适应度清空
%     fitness_sum = zeros(1,population_size);
    % 冒泡排序总适应度
    for i = 1:population_size
        min_index = i;
        for j = i+1:population_size
            if fitness_value(j)<fitness_value(min_index)
                min_index = j;
            end
        end
        % 如果最小值不是当前值，则交换
        if min_index ~=i  
            temp = fitness_value(i);
            fitness_value(i) =  fitness_value(min_index);
            fitness_value(min_index) = temp;
            
            temp_chromosome = population(i,:);
            population(i,:) =  population(min_index,:);
            population(min_index,:) = temp_chromosome;
        end
    end
    for i = 1:population_size
        if i==1
            fitness_sum(i) = fitness_value(i);
        else
            fitness_sum(i) = fitness_sum(i-1) + fitness_value(i);
        end
    end
    % 计算平均适应度
    fitness_avg(generation) = fitness_sum(population_size)/population_size;
    % 保存最佳适应度
    if fitness_sum(population_size) > fitness_best
        fitness_best = fitness_sum(population_size);
        best_generation = generation;
        best_chromosome = population(population_size,:);
    end
    % 记录最佳适应度
    best_fitness(generation) = fitness_best;
    % 逐个选择种群
    new_population=zeros(size(population));
    for i=1:population_size
        % 随机选择种群适应度
        select_fitness_sum = rand * fitness_sum(population_size);
        % 择中法选取最接近随机选择的适应度
        first = 1;
        last = population_size;
        mid = round((first+last)/ 2);
        select_population_id = -1;
        while(first <= last) && (select_population_id == -1)
            if select_fitness_sum > fitness_sum(mid)
                first = mid;
            elseif select_fitness_sum < fitness_sum(mid)
                last = mid;
            else
                select_population_id = mid;
                break;
            end
            mid = round((first+last) / 2);
            if(last-first) == 1
                select_population_id =last;
                break;
            end
        end
        new_population(i,:)=population(select_population_id,:);
    end
    % 更新种群信息
    elitism_population=population(end,:);% 精英种群
    population=new_population;% 选择后的种群
    % 是否保留精英种群
    if elitism
        population(end,:)=elitism_population;
    end
    % 父母交叉染色体(遗传最优值)
    for i = 1:2:population_size
        % 发生交叉
        if rand<cross_rate
            % 随机选择一个交叉位置(可改为多个)
            cross_position = round(rand*chromosome_size);
            % 交叉(互换编码值)
            if (cross_position~=0) && (cross_position~=1)
                temp_cross_chromosome = population(i,cross_position:end);
                population(i,cross_position:end) = population(i+1,cross_position:end);
                population(i+1,cross_position:end) = temp_cross_chromosome;
            end
        end
    end
    % 染色体突变
    for i = 1:population_size
        % 发生染色体突变
        if rand < mutate_rate
            % 随机选择一个交叉位置(可改为多个)
            mutate_position = round(rand*chromosome_size);
            % 变异(编码值取反)
            if (mutate_position~=0)
               population(i,mutate_position)=1-population(i,mutate_position);
            end
        end
    end
end
%% 处理结果
% 解码
max_var = zeros(1,var_num);
for j = 1:var_num
    for k = 1:var_chromosome_size 
        if best_chromosome(1,k+(j-1)*var_chromosome_size) ~=0
            max_var(j) = max_var(j) + 2^(k-1); 
        end
    end
    max_var(j) = lower_bound + max_var(j)*...
                bound_length/(2^var_chromosome_size-1);
end
max_x = max_var(1);
max_y = max_var(2);
max_ans = target(max_x,max_y);
% 打印结果
disp 最优迭代次数
disp(best_generation)
disp 最优个体
disp(best_chromosome)
disp 最优适应度
disp(fitness_best)
disp 最大值位置
fprintf("x=%f\ny=%f\nz=%f\n",max_x,max_y,max_ans);
% 可视化结果
figure('position',[200 200 1280 480])
subplot(1,2,1);hold on;
[x,y]=meshgrid(-4:0.1:4);
z=target(x,y);
mesh(x,y,z);
title('Target Function')
xlabel('x');ylabel('y');zlabel('z');
view(-37.5,30)
pause(0.1)
scatter3(max_x,max_y,max_ans,'r*');
subplot(1,2,2);hold on;
x = 1:generation_size;
y1 = fitness_avg;
y2 = best_fitness;
try
    assert(length(x)==length(y1))
    assert(length(x)==length(y2))
    yyaxis right
    plot(x,y1,'-','LineWidth',1.2)
    yyaxis left
    plot(x,y2,'-','LineWidth',1.2)
    xlabel('x');ylabel('y');
    legend('fitness avg','best fitness')
    title('迭代数据')
catch
    disp("plot failed")
    disp(size(x))
    disp(size(y1))
    disp(size(y2))
end

% function y=optimization_target(x)
%     y=x+10*sin(5*x)+7*cos(4*x);
% end